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Abstract 

Background: MicroRNAs (miRNAs) play important roles in regulating the expression of protein-coding genes by 
directing the degradation and/or repression of the translation of gene transcripts. Growing evidence shows that 
miRNAs are indispensable player in organismal development with its regulatory role in the growth and 
differentiation of cell lineages. However, the roles of miRNA-mediated regulation in environmental adaptation of 
organisms are largely unknown. To examine this potential regulatory capability, we characterized microRNAomes 
from the brain of zebrafish raised under normal (28°C) and cold-acclimated (10°C, 10 days) conditions using Solexa 
sequencing. We then examined the expression pattern of the protein-coding genes under these two conditions 
with Affymetrix Zebrafish Genome Array profiling. The potential roles of the microRNAome in the transcriptomic 
cold regulation in the zebrafish brain were investigated by various statistical analyses. 

Results: Among the total 214 unique, mature zebrafish miRNAs deposited on the miRBase website (release 16), 
175 were recovered in this study. In addition, we identified 399 novel, mature miRNAs using multiple miRNA 
prediction methods. We defined a set of 25 miRNAs differentially expressed under the cold and normal conditions 
and predicted the molecular functions and biological processes that they involve through Gene Ontology (GO) 
annotation of their target genes. On the other hand, microarray analysis showed that genes related to mRNA 
processing and response to stress were overrepresented among the up-regulated genes in cold-stress, but are not 
directly corresponding to any of the GO molecular functions and biological processes predicted from the 
differential miRNAs. Using several statistical models including a novel, network-based approach, we found that 
miRNAs identified in this study, either individually or together, and either directly or indirectly (i.e., mediated by 
transcription factors), only make minor contribution to the change in gene expression patterns under the low- 
temperature condition. 

Conclusions: Our results suggest that the cold-stress response of mRNA expression may be governed mainly 
through regulatory modes other than miRNA-mediated regulation. MiRNAs in animal brains might act more as 
developmental regulators than thermal adaptability regulators. 



Background 

MicroRNAs (miRNAs) are a class of -22 nt, non-coding 
RNAs present in plants and animals that play important 
roles in regulating gene expression by directing transla- 
tional repression and/or mRNA destabilization [1-3]. 
Recognition and interaction between miRNA and 
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mRNA take effect only if the 5' end seed (6-8 nt) of the 
miRNA is complementary to the 3' UTR site of the 
mRNA. A miRNA typically targets multiple, even hun- 
dreds of, genes [4], which makes miRNAs potentially 
global regulators responsible for spatio-temporal gene 
expression patterns. 

Like let-7 [5,6], some miRNAs are highly conserved 
throughout evolution and participate in pivotal biologi- 
cal processes [7]. Lineage- or species-specific miRNAs 
may contribute to the uniqueness of organisms [8,9]. 
With recent advancements in miRNA microarrays and 
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next-generation sequencing, miRNA expression profiles 
of many species are being generated at a rapid pace. 
These profiles are invaluable resources for better under- 
standing gene expression dynamics [10-16]. Some miR- 
NAs have been well characterized and have been shown 
to have important functions in normal development 
[17-21]. The deficient expression of some miRNAs also 
presents severe phenotypic abnormalities [22,23]. The 
role of miRNAs in responses to cellular and environ- 
mental stress has been discovered in a variety of organ- 
isms [24-28]. A few of these miRNAs have been 
characterized as key regulators of the stress response 
[25]. For example, miR-8 family miRNAs were shown to 
be indispensable in the response of zebrafish embryos to 
osmotic stress [29]. Despite this progress, most studies 
to date have only examined miRNAs' functions at indi- 
vidual level. The global contributions of microRNAomes 
(miRNAomes) to environmental adaptations remain lar- 
gely unexplored. 

The zebrafish, an important model organism, is an 
eurythermal fish with temperature tolerance ranging 
from approximately 7°C to 40°C [30], making it a useful 
experimental subject for exploring the molecular 
mechanisms underpinning thermal responses. Almost all 
of the previous work on the molecular underpinnings of 
thermal responses has focused on the fluctuated expres- 
sion of protein-coding genes in response to warm or 
cold acclimation [31,32]. Many miRNAs have been dis- 
covered in zebrafish [33,34], and their expression pro- 
files [35] and functions in developmental regulation 
have been studied [36]. Our goal was to determine how 
the zebrafish miRNAome responds to thermal acclima- 
tion and to assess the roles of miRNAs in thermally 
challenged transcriptomic plasticity. Toward this end, 
we systematically characterized the miRNAomes of zeb- 
rafish brains under normal and cold-acclimated condi- 
tions using Solexa sequencing. We also identified 
genome-wide gene expression patterns under each of 
these conditions using microarray profiling. We then 
examined the potential links between the changes in 
miRNA levels and the transcriptomic shifts and found 
that miRNAs may play minor roles in transcriptomic 
plasticity during cold acclimation in the zebrafish brain. 

Results 

Deep sequencing of the zebrafish brain small RNA and 
miRNA prediction 

Total RNAs from brain of cold-treated (10°C, 10 days) 
and normal control (28°C) zebrafish were extracted and 
electrophoresed to enrich the 18-30 nt small RNAs and 
then subjected to Solexa sequencing. The numbers of 
high-quality reads were 9,615,499 and 14,399,653 for the 
normal control (NC) and cold-acclimated (WC) RNA 
libraries, respectively. The length distributions of the 



sequence reads from the two samples both peaked at 22 
nt, conforming to size distribution patterns characteris- 
tic of miRNA libraries (Figure 1). After filtering short 
(<18 nt) and adaptor-contaminated reads, we obtained 
6,473,829 (NC) and 10,799,032 (WC) redundant reads 
that corresponded to 349,107 (NC) and 350,474 (WC) 
unique reads. 

We ran the SOAP2 program [37] to map the unique 
reads to the zebrafish genome (Danio_rerio.Zv8.56.dna. 
toplevel.fa.gz). Stem-loop structures were predicted from 
the putative miRNA precursor sequences using MIREAP 
(https://sourceforge.net/projects/mireap/). Overall, we 
identified 602 stem-loop precursors that coded for 700 
miRNAs. Among them, 175 were previously known 
miRNAs [38], constituting over 80% of the unique miR- 
NAs (214) deposited in the miRBase database (release 
16). 525 are newly identified miRNAs (Additional file 1). 
To cross-validate these newly identified miRNAs, we 
adopted another miRNA prediction software, miRDeep 
[39] by invoking the quantifier.pl script (using default 
parameters except that -g was set to 0) in the software 
package. Results showed that the miRNA prediction 
outcomes are almost the same, independent of the pre- 
diction methods used and the miRNA expression levels 
are also highly consistent between the two methods 
(Additional files 2, 3 and 4). To further evaluate the 
accuracy of miRNA prediction, we examined the pre- 
dicted miRNAs with a higher stringency by adding 
another two criteria: 1) they cannot be exon-derived (a 
recent study by Stark et al [40] suggested that miRNAs 
lie almost exclusively in the introns and intergenic 
regions of fruit flies); 2) their precursors need to be 
recognized as genuine stem-loop structures by the 
MiPred [41] software. Consequently, the number of 
newly identified miRNAs in our dataset was reduced 
from 525 to 399 (See Additional file 1). However, 
MiPred does not ensure 100% accuracy in determining 
the authenticity of miRNAs. For instance, 5 of the 175 
well-known zebrafish miRNAs were incorrectly assigned 
in this prediction. In addition, "exon-derived" miRNAs 
are possible to exist if the relevant exons do not code 
for proteins or are transcribed under alternative promo- 
ters or from alternative strand. Therefore, the number 
of authentic miRNAs in the zebrafish brain is very likely 
higher than 574 (i.e., 175+399). To verify the authenti- 
city of the newly identified miRNAs, we selected 13 
miRNAs expressed at various levels to perform RT-PCR. 
All of them can be successfully amplified (Additional 
file 5). Cloning and sequencing of the PCR products 
confirmed that they are in correct forms. These results 
showed that the newly identified miRNAs are very likely 
true miRNAs. Remarkably, the newly identified miRNAs 
are expressed at substantially lower levels than pre- 
viously known miRNAs (Figure 2). This may partially 
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Figure 1 Length distribution of the high-quality reads in the 10-30 nt range. 



explain why the new miRNAs were missed by traditional 
cloning and sequencing methods. 

Solexa sequencing are known for its short reads, and 
errors can be pervasive in short-read sequencing [42]. 
To address this concern, we clustered the raw reads 
using the FreClu software [43], which takes the read 
count into consideration. The results showed that this 
clustering process can improve the accuracy of short 
reads mapping (Table 1). However, it did not improve 
the data quality of the predicted miRNAs (Additional 
file 6). Hence, we performed all of the follow-up ana- 
lyses of the miRNAs identified using the original reads. 

The miRNA expression profile in cold acclimation 

Based on normalized read counts of each miRNA in the 
control and cold libraries, we deduced the expression 
changes of the brain miRNAome in response to cold 
acclimation. Because low-expressed sequences have 
more noise than high-expressed sequences, Koh et al. 
proposed an adaptive thresholding method using Kolo- 
mogorov-Smirnov statistics to determine the threshold 
value for identifying significant miRNAs [44]. We 
observed a similar pattern in the present study. The less 
abundant miRNAs showed a greater divergence than the 



more abundant miRNAs (Figure 3). We hypothesized 
that the low-expressed miRNAs showing greater diver- 
gence likely have more noise than the high-expressed 
ones and that the high-expressed miRNAs have a high 
probability of playing an important biological role. We 
borrowed the concept of Koh et al. to define miRNAs 
differentially expressed under the normal and cold con- 
ditions as the ones that have a greater than 2-fold 
change between the two libraries and a mean expression 
level greater than 64. As a result, we found 11 up-regu- 
lated and 15 down-regulated miRNAs in the cold 
library; a majority of them are known miRNAs (Table 
2). 

To reveal which biological processes might be regu- 
lated by the changed miRNAome, we predicted the tar- 
get genes of the differential miRNAs using PITA [45]. 
We then projected the functions of these differential 
miRNAs onto the overrepresented Gene Ontology (GO) 
terms of their targets using the GO:TermFinder package 
[46] with a corrected p value of <0.05 as a cutoff. In 
total, we found 4 (dre-mir-30d-5p, dre-mir-182-5p, dre- 
mir-458-3p, and miR_8) and 5 miRNAs (miR_46, dre- 
mir-133a-2-3p, dre-mir-146a-5p, dre-mir-183-5p, and 
miR_8) with overrepresented biological process and 
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Figure 2 Violin plot showing the expression levels of novel-identified and known miRNAs. The violin plot combines box plot and kernel 
density trace to describe the distribution pattern of a vector of data. The four violin plots were drawn according to the expression levels of 525 
newly identified and 175 known miRNAs under the normal and cold-acclimated condition, respectively. Note that if the expression level of a 
miRNA was 0 under certain condition, we assigned it to 1 in order to make the log-transfer of the raw data. 



molecular function GO terms, respectively (Table 3). In 
addition, we assembled all target genes of the up-regu- 
lated and down-regulated miRNAs and examined 
whether any overrepresented functions would emerge. 
Our results showed that genes related to "purine 
nucleotide binding" and "nucleotide binding" were 
enriched in the targets of up-regulated miRNAs. 



Gene expression profiling and differentially expressed 
genes 

Gene transcription changes between the cold-accli- 
mated and normal samples were obtained by microar- 
ray-based hybridization using three sets of RNA 
samples prepared from fishes of three acclimation 
experiments. The gene expression level was measured 



Table 1 SOAP mapping results of reads with or without clustering. 



Before Clustering After Clustering 







unique 


% 


redundant 


% 


unique 


% 


redundant 


% 


NC 


# reads 


349107 




6473829 




255724 




6473829 






# matched 


102810 


29.45 


5016748 


77.49 


99235 


38.81 


5206070 


80.42 


WC 


# reads 


350474 




10799032 




287325 




10799032 






# matched 


153533 


43.81 


8680394 


80.38 


151896 


52.86 


8861106 


82.05 



WC: cold-acclimated condition; NC: normal condition. 



Yang et al. BMC Genomics 201 1, 12:605 
http://www.biomedcentral.eom/1 471-21 64/1 2/605 



Page 5 of 1 7 




Expression level of miRNAs in cold library (log scale) 



Figure 3 Scatter plot showing the miRNA expression levels of the zebrafish cold and normal brain libraries. The plot indicates much 
more expression noise of lowly expressed miRNAs than that of highly expressed ones. Hence, to define biologically significant miRNAs, we 
screened the differential genes from those with mean expression level greater than 64. 



as the intensity of probe sets on an Affymetrix micro- 
array. After filtering the probe sets and retrieving the 
map relationship between the probe sets and the 
genes, we obtained gene expression profile recordings 
of 7356 protein-coding genes. The differential genes 
were determined using a two-class, unpaired compari- 
son in the SAM (version 3.0) program [47]. Using a 
5% false discovery rate (FDR) cut-off, we found that 
311 genes were down-regulated and 315 genes were 
up-regulated in the cold-acclimated samples (Addi- 
tional file 7). There were 127 and 85 genes with more 
than 2-fold changes in mRNA abundance that were 
down- and up-regulated, respectively. 

We then investigated the functional enrichment of 
these differential genes using the GO:TermFinder 



package [46]. As shown in Table 4, under cold stress, 
the genes with a functional annotation related to 
"mRNA processing" were enriched in the up-regulated 
differential genes, and genes associated with "cellular 
macromolecular complex subunit organization" were 
overrepresented in the down-regulated genes. An inter- 
esting observation here is that the identified molecular 
functions and biological processes changes in the cold- 
acclimated zebrafish brain bear no apparent overlap 
with those predicted from the potential targets of the 
differential miRNAs (see Table 3), suggesting the role of 
miRNAs might be minor in the transcriptomic changes 
in cold stress. We thus examined the possible regulation 
roles of the microRNAome in cold adaptation by other 
statistical approaches. 
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Table 2 Differentially expressed miRIMAs of zebrafish brain under cold-acclimated and normal conditions. 



miRNA name 


Genome coordinate 


miRNA sequence 


# wc 


# NC 


Fold A 


d re-m i r-9-5-3 p 


5:49 ] 44568:49 1 44646:- 


UAAAGCUAGAUAACCGAAAGU 


37364 


2559 


7.377 (|) 


dre-mir-30d-5p 


1 6:238601 43:23860227:+ 


UGUAAACAUCCCCGACUGGAAGCU 


19903 


3770 


2.667 (|) 


miR 49 


4:52193679:52193757:+ 


UUCACUGUGGCGGAAAUGACC 


146 


31 


2.380 (|) 


miR 46 


21 :25385786:25385876:+ 


AAGAGAAGAGUGAGCGAGUGA 


149 


32 


2.352 (f) 


dre-mir-30e-2-5p 


1 3:28550257:28550337:+ 


UGUAAACAUCCUUGACUGGAAGC 


2856 


617 


2.339 (t) 


d re-m i r-24-4-3 p 


22:5036700:5036779:+ 


UGGCUCAGUUCAGCAGGAACAG 


652 


142 


2.320 (|) 


dre-mir-21 91 -3p 


20:5 1 1 08008:5 1 1 08091 :+ 


UCACACCUACAAUCCCUGGCA 


435 


99 


2.221 (t) 


dre-mir-1 46a-5p 


13:1 1738901:1 1738980:+ 


UGAGAACUGAAUUCCAUAGAUGG 


42728 


10063 


2.145 (t) 


miR 48 


24:36888097:368881 73:+ 


UGAGGAGUUUAGAGCAAGUAA 


129 


31 


2.102 (t) 


miR 23 


5:8733083:87331 58:- 


AGCUGGUGUCCUGCAGAGUUU 


345 


85 


2.051 (|) 


d re-m ir-45 5-5 p 


5:66430757:66430834:- 


UAUGUGCCCUUGGACUACAUC 


168 


42 


2.021 (t) 


dre-mir-30b-3p 


1 6:23860432:238605 1 0:+ 


CUGGGCGGAGGGUGUUUGCU 


456 


462 


2.006 (j) 


dre-mir-1 45-5p 


1 4:4045 1 735:4045 1816:+ 


GUCCAGUUUUCCCAGGAAUCCC 


104 


106 


2.017 (I) 


dre-mir-1 29-1 -3p 


4:1 8861 279:1 8861 364:- 


AAGCCCUUACCCCAAAAAGUAU 


643 


697 


2.145 (j) 


dre-mir-1 33a-2-3p 


23:5570884:5570961:- 


UUGGUCCCCUUCAACCAGCU 


117 


130 


2.199 (I) 


miR 13 


5:6754333:6754409:+ 


UGUGGUGAUUAGAACAGGCUGA 


159 


181 


2.253 (j) 


miR_3 


Zv8_NA1 0078:5 11 97:5 1280:- 


AUGACUCAAACCCGAGGACUU 


1279 


1618 


2.504 (|) 


dre-mir-92a-1-3p 


8:45338553:45338629:+ 


UAUUGCACUUGUCCCGGCCU 


9007 


11736 


2.579 (|) 


dre-mir-458-3p 


14:26650559:26650634:- 


AUAGCUCUUUGAAUGGUACU 


17139 


22528 


2.601 (|) 


dre-mir-1 29-2-3p 


7:50935695:50935777:- 


AAGCCCUUACCCCAAAAAGCAU 


7747 


10824 


2.766 (|) 


dre-mir-1 52-3p 


3:21062159:21062235:+ 


UCAGUGCAUGACAGAACUUUG 


2271 


3467 


3.021 (|) 


dre-mir-1 83-5p 


4:14884980:14885059:+ 


UAUGGCACUGGUAGAAUUCACU 


182 


303 


3.294 (|) 


miR_5 


Zv8_NA3780:1 1797:1 1893:+ 


GCAUUGGUGGUUCAGUGGGA 


698 


1247 


3.536 (|) 


dre-mir-1 82-5p 


4:14886067:14886150:+ 


UUUGGCAAUGGUAGAACUCACA 


151 


365 


4.783 (|) 


miR_8 


2:35918353:35918422:+ 


AGCGGAUCAACUGAUGCGAC 


98 


394 


7.956 (|) 



#WC: the number of miRNAs in the cold library; #NC: the number of miRNAs in the normal library; Fold A : normalized fold change of miRNAs in the normal over 
the cold library. |(|) denotes the up(down)-regulation. 



Verification of the expression pattern of the protein- 
coding genes and miRNAs by quantitative RT-PCR (q RT- 
PCR) 

The reliability of statistical analyses depends on the 
quality of the microarray and solexa profiling data. To 
evaluate the accuracy of the microarray screening and 
the miRNA profiling results, genes and miRNAs with 
change direction either down-regulated (plplb, h3f3c 
and dre-mir-24-4-3p) or up-regulated (ernx2, ikzf5, 
dazap, mmpl3, zgc:112425, taldol, and dre-mir-1 52-3p) 
were selected for quantitative RT-PCR verification. Q- 



PCR showed that the direction of expressional change 
for each selected gene correlated perfectly with the 
results of the microarray screening or solexa sequencing 
(Figure 4). Although the changed folds of expression in 
a few genes {mmpl3a, plplb) and the dre-mir-152-3p 
showed some discrepancy between the q-PCR results 
and the '-omic' (microarray and solexa sequencing) pro- 
filing data, the other genes/miRNAs are in good agree- 
ment. Since we evaluate only the direction of change of 
the genes and miRNAs when delineating the potential 
relationship between miRNA-mediated regulation and 



Table 3 Differentially expressed miRNAs and their functions. 


miRNA name 


FoldA 


Molecular Function 


Biological process 


dre-mir-30d-5p 


2.667 (t) 


ND 


Phosphoinositide metabolic process 


miR_46 


2.352 (t) 


ribonucleotide binding, ATP binding 


ND 


dre-mir-1 46a-5p 


2.145 (T) 


protein tyrosine/serine/threonine phosphatase activity 


ND 


dre-mir-1 82-5p 


4.783 (|) 


ND 


glucose metabolic process 


dre-mir-458-3p 


2.601 (|) 


ND 


regulation of developmental process 


miR_8 


7.956 (|) 


aminoacyl-tRNA ligase activity 


tRNA aminoacylation for protein translation 


dre-mir-1 33a-2-3p 


2.199 (|) 


receptor signaling protein activity 


ND 


dre-mir-1 83-5p 


3.294 (|) 


cytoskeletal protein binding 


ND 



ND: Not determined. Fold A : fold change of miRNAs in the cold over the normal library. denotes the up(down)-regulation. 
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Table 4 Functional enrichment of differentially expressed mRNAs. 


Category 


# genes 


Molecular Function 


Biological process 


Up-regulated 


315 


nucleotide binding 


mRNA processing 


Down-regulated 


311 


ND 


cellular macromolecular complex subunit organization 


Up-regulated & >2-fold change 


85 


ND 


response to stress 








response to temperature stimulus 


Down-regulated & >2-fold change 


127 


dipeptidase activity 


cellular macromolecular complex assembly 






oxidoreductase activity 


nucleosome assembly 



# genes: the number of genes in the category; ND: not determined. 



transcriptomic shifts (see later sections), the gene and 
miRNA expression profiles obtained from the '-omic' 
profiling approaches are generally useful for statistical 
analyses. 

Statistical tests for the potential correlation between the 
expression levels of miRNAs and protein-coding genes 

Mounting evidence has shown that miRNAs can orches- 
trate the expression of mRNAs and that the well- 
balanced interplay between the miRNAome and protein- 
coding gene expression plays a crucial role in an organ- 
ism's development [48-50]. To investigate the influence 
of miRNAs on the global expression of protein-coding 
genes in cold stress, we adopted four different 
approaches. First, we estimated the effect of miRNAs on 
differentially regulated genes using Fishers exact test. 
Using PITA target prediction methods, we obtained 



18665 pairs of miRNA-target interactions in total. Our 
results showed only one miRNA (dre-mir-34-3p) with 
target genes overrepresented in the down-regulated gene 
set. No miRNAs showed targets enrichment or depletion 
in the up-regulated gene set (FDR < 5%). The miRNA 
dre-mir-34-3p was expressed at very low levels (29 in 
the cold library and 20 in the control library) in the 
brain and was not significantly different between the 
two conditions. We speculate that it is not probable that 
dre-mir-34-3p plays an important role in the changed 
gene expression patterns observed under cold stress. To 
address concerns over potential target prediction bias by 
PITA, we performed target prediction by miRanda soft- 
ware [51]. Under this methodology, we obtained 56670 
pairs of miRNA-target interactions, which is about 3- 
fold of the number (18665) obtained by PITA prediction 
method. We then examined the composition of targets 
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Figure 4 Fold change of expression level for selected genes and miRNAs measured by q RT-PCR. Quantitative RT-PCR shows that the 
direction of expressional change under the cold condition of the selected genes and miRNAs correlated nicely with the results obtained from 
the "-omics" data (Additional file 7 and Table 2). 
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predicted by the two methods for each miRNA. Our 
results showed that 615 miRNAs shared at least one 
predicted target between PITA and miRanda; however, 
as a whole, the overlap is relatively low (Additional file 

8, Table S3-S6). For example, among a group of 370 
miRNAs for which each has more than ten targets by 
whatever target prediction methods, the median value of 
the proportion of overlap is only 0.396 (Additional file 

9, Figure S3). In average, each miRNA targets -80 pro- 
tein-coding genes and each gene is likely to be regulated 
by 7-8 miRNAs in return. Based on these interaction 
data, our statistical analyses recaptured the fore-men- 
tioned miRNA-mRNA interacting pattern, i.e., no miR- 
NAs targeted over- or under-represented genes in either 
up- or down- regulated gene set. 

Recently, Chen et al. [52] proposed a novel concept 
for the regulatory effect score (RE score) to measure the 
inhibition capability of a miRNA on its targets. In prin- 
ciple, the RE score is calculated as the average difference 
in expression levels between a target and a non-target 
group. By calculating the RE scores on the PITA and 
miRanda methodology, respectively, we checked whether 
some of the miRNAs would changed regulatory effects 
on their targets between the two conditions. Consistent 
with previous results, we again found that no miRNAs 
exhibited a significant change in their regulatory effect 
after correcting the raw t-test p-values by the multiple 
hypothesis test (FDR < 0.05) whatever the target predic- 
tion method was used. Taken together, our results indi- 
cate that individual miRNAs are unlikely to be major 
factors leading to the changed transcriptomic phenotype 
of zebrafish brain under cold stress. 

To further reduce the concerns over target-prediction 
bias induced errors, we applied the Sylamer algorithm 
[53] to re-examine the potential relationship between 
the miRNAome and transcriptome under cold acclima- 
tion. We sorted the protein-coding genes from the most 
significantly up-regulated to the most down-regulated. 
We then calculated the distribution of the enrichment/ 
depletion of words in the 3' UTR sequences that could 
be complementary to the seed sequences of miRNAs. 
Since positions 1-7 and 2-8 of the 5' end of a mature 
miRNA are the most important sequences for target 
recognition [54], we took the 7 nt region as a miRNA 
seed in this study and examined whether some over- or 
under-represented words corresponded to certain 
miRNAs. 

As shown in Figure 5, there are no obvious peaks in 
the region containing the most down-regulated genes, 
suggesting that the analysis failed to identify any miR- 
NAs being the major and direct factors contribute to 
the decreased expression levels of the down-regulated 
gene sets; on the other hand, few than ten words seem 
to be overrepresented in the 3' UTR regions of the up- 



regulated genes. The three most overrepresented words 
are TAATAGA (P = 3.86 x 10" 10 ), GCGAATT (P = 1.58 
x 10" 9 ) and GTTCACG (P = 4.58 x 10" 9 ). However, 
these words are not complementary to any seed region 
of the profiled zebrafish miRNAs. We then examined 
the composition of the top 10 most enriched words and 
found only one, GCTAGTT (P = 5.53 x 10" 8 ), that was 
complementary to the seed region of the newly identi- 
fied miR_381 (AACTAGCAGCTGTTGACATCCA). 
However, this miRNA was hardly expressed in the brain 
(1 in the control library and 4 in the cold library), sug- 
gesting that it was unlikely to be the major factor 
responsible for the changed gene expression patterns 
observed. 

To systematically investigate the extent to which miR- 
NAome-mediated regulation could have influenced the 
transcriptome under cold-acclimation, we designed a 
novel, network-based methodology to parse the subtle 
regulatory effects of a group of miRNAs on a group of 
protein-coding genes at the system level. Based on the 
predicted miRNA-target relationship predicted from the 
PITA method, we first built a miRNA-mRNA network 
of 32 differential miRNAs and 626 mRNAs (311 down- 
regulated and 325 up-regulated). The miRNA-target 
relationship showed that the majority of the differential 
gene transcripts do not possess any potential targeting 
sites in the 3' UTRs that correspond to the differentiated 
miRNAs. In addition, approximately half of the differen- 
tial miRNAs do not target any genes in the differen- 
tiated gene list. The resultant miRNA-mRNA network 
involved only 97 nodes (82 mRNAs and 15 miRNAs) 
and 97 edges when the isolated nodes were separated. 
As seen in Figure 6, the miRNA-mRNA interaction net- 
works are composed of six unconnected components. 
One component contains 60 of the 97 nodes, whereas 
the others are formed by only 1 or 2 miRNAs situated 
in the network center. The node with the highest num- 
ber of edges is the up-regulated miRNA miR_46 that 
was predicted to target 18 (11 up-regulated and 7 
down-regulated) gene transcripts. Among the 82 
miRNA-targeted differentially expressed genes, over 
four-fifths (67/82) are predicted to be targeted by only 1 
miRNA and one-fifth by 2 or more. We then computed 
the "regulatory density" (RD) of the network by formula 
(1) (see Methods). We developed this formula and 
define it as a proxy for the regulatory potential of a 
group of miRNAs upon gene transcripts by calculating 
the proportion of coherent over incoherent regulatory 
links within a miRNA-mRNA network. Our analysis 
showed that the observed RD was 0.1134. However, the 
values of the 95th- and 99th-quantiles were 0.1546 and 
0.1959, respectively, for the 1000 cohorts of random net- 
works that were constructed by shuffling the original 
miRNA-mRNA network while maintaining the degree 
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distribution. We further analyzed the topologic property 
of a similar network derived by using the miRanda pre- 
diction data. The produced miRNA-mRNA network 
consisted of 203 edges that link 25 and 161 differential 
miRNAs and protein-coding genes, respectively. As seen 
in Additional file 10, a majority of the mRNAs are tar- 
geted by a few miRNAs; and only a few miRNAs have 
>5 target genes. Calculation showed that the RD score 
of this network is 0.0148, again not significantly greater 
than random expectation (empirical P-value: 0.269). 
Taken together, results from all the above statistical 
analyses suggested that the miRNAome is unlikely to be 
the major regulatory component responsible for the 
changed expression profiles of the protein-coding genes 
under cold stress. 

Mathematical models [55] and statistical associations 
[56,57] have shown that miRNAs frequently intertwine 
with transcription factors (TFs) in regulatory circuits 
and play essential roles in signal transduction. 
Recently, Tu et al. [58] proposed a two-layer network 



wherein miRNA-initiated regulation was mediated and 
propagated through the miRNAs' targeted transcription 
factors. We investigated whether the miRNAs had 
indirectly, through their transcription factors, led to 
the observed transcriptomic plasticity. We first filtered 
the transcription factors from the arrayed gene list 
according to the GO annotation. We treated the genes 
with any one of the GO terms, "DNA binding", "tran- 
scription factor activity", "transcription activator activ- 
ity" or "transcription repressor activity", as putative 
transcription factors. Accordingly, we predicted 796 
TF-coding genes among the 7356 protein-coding 
genes. Notably, only 46 of these TFs showed differen- 
tial expression in cold acclimation, and they were sig- 
nificantly underrepresented in the differential 
expressed gene list (46 out of 626 differential genes 
versus 796 TFs in the 7356 arrayed genes, Fisher's 
exact test: P < 0.01). Among the differential TF genes, 
only 6 were predicted to be the targets of 5 differential 
miRNAs (the pair in Figure 6). Of the 6 TFs, 3 



Yang et al. BMC Genomics 201 1, 12:605 
http://www.biomedcentral.eom/1 471 -21 64/1 2/605 



Page 1 0 of 1 7 




Figure 6 Graphic illustration of the regulatory relationship between the differential miRNAs and mRNAs predicted by PITA. Triangles, 
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miRNAs or mRNAs, respectively. The edges in dark and yellow denote coherent and incoherent regulation, respectively. The size of a node is 
proportional to the number of its partners. Of nodes with degree > 5, the associated name is shown. 
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showed coherent regulatory direction with their 
respective miRNA partners (miR_13^ emx2; dre-mir- 
183-5p-H'fe/5; miR_46, miR_48->/z3/3c. see Figure 6). 
Given that the gene emx2 represented slightly up-regu- 
latoin (fold-change: 1.214; q-value: 4.376%) under cold- 
acclimated condition, we speculate that it was not a 
driven gene responsible for the changed mRNA 
expression pattern. There is no evidence that the ikzfS 
(fold-change: 1.994; q-value: 0) and h3f3c (fold-change: 



2.245; q-value: 2.144%) are related to stress response 
when we examined their annotations in multiple data- 
bases and in the literature. Additionally, almost no dif- 
ferential TF-coding genes were targeted by differential 
miRNAs. Taken together, we conclude that the mRNA 
phenotypic plasticity seen under the cold environment 
was not dominantly modulated by miRNAs, even if 
these factors were directly or indirectly mediated by 
transcriptional factors. 
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Discussion 

The repertoires of miRNAs have been constantly revised 
and updated with the advance of in silico prediction 
methods and more intensive cloning and sequencing 
efforts in broader tissues and developmental stages. The 
number of unique, mature human miRNAs is now up 
to 1223 in the miRBase database (version 16). In con- 
trast, there are only 214 mature miRNAs now registered 
for zebrafish. This relatively low number could be the 
result of genomic between-species differences or incom- 
plete miRNA identification in zebrafish. By Solexa 
sequencing, we comprehensively profiled the miRNA 
repertoire in the zebrafish brain. The number of pre- 
dicted miRNAs in this tissue reached 700 (574 in high- 
stringency criteria) in this study and greatly expanded 
the miRNAome in zebrafish. Of note is that the expres- 
sion levels of the newly identified miRNAs in the zebra- 
fish brain are mostly low. The authenticity of these 
newly identified miRNAs remains to be further investi- 
gated to determine their biological functions. However, 
a population analysis on 7 randomly selected novel miR- 
NAs currently carried out in the laboratory showed that 
almost half of them are subjected to purified selection 
suggesting that many newly identified miRNAs may 
have high potentials to be functional in vivo despite 
their low expression levels. It is possible that some of 
these newly identified, low-expressed miRNAs could be 
expressed at higher levels in tissues other than the brain 
and play roles in gene regulation. The miRNAome of 
the zebrafish is much larger than those of the lower 
chordates such as the amphioxus species, which had 
been determined by small RNA library sequencing [59] 
and solexa sequencing [13]. This suggests the expansion 
of the miRNA regulome in the evolution of vertebrates. 

To study the function of the miRNAome in environ- 
mental adaptation, we designed the cold-acclimation 
experiments in this study. We profiled both the miR- 
NAome and the transcriptome of the brain under low- 
and normal-temperature conditions. We adopted four 
approaches to examine how miRNAs, either individually 
or globally, could have influenced the transcriptomic 
plasticity under cold conditions. Using statistical meth- 
ods we first examined the potential effects of individual 
miRNAs on gene expression patterns. We found one or 
two miRNAs exhibiting regulation signatures on mRNA 
expression based on enrichment/depletion analysis of 
target genes among the differential gene set. However, 
we argue that these miRNAs are unlikely to have been 
responsible for the different expression patterns after 
further investigation of their expression levels. Our 
results show that at the individual level, the miRNAs 
were not a main factor shaping the mRNA expression 
patterns observed under the cold-acclimated conditions. 
We then attempted to clarify whether the cold-stress 



response of mRNA expression is governed by miRNAs 
indirectly and is mediated by TFs. TF-mediated tran- 
scriptional regulation is thought to be stronger than the 
post-transcriptional regulation usually mediated by miR- 
NAs [60]. For this purpose, we assigned a set of TFs 
according the GO annotation and found that TFs were 
significantly underrepresented in the differential protein- 
coding genes; there were only 6 TFs targeted by differ- 
ential miRNAs. Again, our results suggest that there is 
hardly any influence from miRNAs on the changed 
mRNA expression patterns. A recent study by Nunez- 
lglesias et al [50] suggested that a subtle regulation rela- 
tionship between miRNAs and their targets can be 
observed in the global miRNAs-mRNA network using 
more sensitive experimental designs and a systems biol- 
ogy approach. Aiming to discern the subtle regulatory 
effects, we employed a novel network-based analysis. 
We found that the miRNA-mRNA network (formed by 
the regulatory relationship between differential miRNAs 
and differential mRNAs) does not possess a significantly 
higher percentage of coherent edges than the random 
expectation. Furthermore, this regulatory network con- 
sists of only 3 pairs of coherent edges between miRNAs 
and TF genes. Therefore, our data also rule out the pos- 
sibility that miRNAs as major factors indirectly modu- 
late mRNA expression in the brain of cold-acclimated 
zebrafish. One might argue that many TFs are often 
expressed in low abundance and that some of them 
might go beyond the detective sensitivity of the current 
Affymetrix array screening methodology, causing the 
observed underrepresentation of the TFs in the differen- 
tial gene set. This situation might eventually lead to 
incomplete evaluation of the regulatory effects of miR- 
NAs on the transcriptome. To evaluate this possibility, 
we calculated the mean expression level of the 796 pre- 
dicted TFs and the non-TF genes in the array. As seen 
from figure 7, though the TFs showed slightly lower 
average expression levels than the non-TF genes, the 
distribution pattern of expression intensity for TFs as 
compared with non-TFs within the differential genes is 
very similar as that within the whole set of the arrayed 
genes, indicating that the relatively low expression of 
TFs did not hinder the screening of the differential TF 
genes in this study. Evidence from various developmen- 
tal biology studies has been continuing emerged favor- 
ing the hypothesis that miRNAs are capable of 
imparting robustness on gene regulation networks that 
involve several classes of recurrent circuits composed of 
miRNAs and their corresponding targets [61]. For exam- 
ple, a study by Hilgers et al [62] revealed that miR- 
263a/b plays an important role in protecting Drosophila 
mechanosensory bristles from apoptosis and ensures the 
developmental robustness of this organ. In our study, 
we found miRNAs are minor factors in transcriptome 
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genes is very similar as that within the whole set of the arrayed genes, indicating that the relatively low expression of TFs did not influence the 
screening of the differential TF genes in this study. DEG: differential genes; non-DEG: non-differential genes. 



changes in the cold acclimated zebrafish brain. This dif- 
ference may reflect that miRNAs in animal brains act 
more as developmental regulators than adaptability 
regulators. 

Thermally invoked responses can be generally classi- 
fied as cellular stress responses, homeostatic responses 
or body remodeling-related responses that correspond 
to short-term, mid-term or long-term acclimation pro- 
cesses, respectively [63-65]. Our samples were taken 
after 10-day cold acclimation, which can be classified as 
mid-term acclimation. We reason that the current 



results reflect only the gene regulatory patterns of the 
zebrafish brain at this temporal scale of acclimation. 
The involvement of miRNAs as a pivotal regulator in an 
earlier or later cold adaptation phase remains to be 
further investigated. 

Persistent efforts have identified and characterized 
some cold stress-induced protein-coding genes [66,67]; 
among them, rbrn3, which encodes a glycine-rich RNA 
binding protein, has been suggested to regulate global 
protein synthesis by affecting miRNA expression levels 
[68]. In the present study, we also observed two 
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homologous genes in zebrafish, cirbp (fold-change: 
2.244; q-value: 0) and zgc:112425 (fold-change: 1.961; q- 
value: 1.598%), with significantly up-regulated expression 
levels under cold-acclimated condition. 

Lopez-Maury et al [69] have pointed out that post-tran- 
scriptional regulatory mechanisms, compared to transcrip- 
tional mechanisms, are largely unknown under stress 
conditions. To our knowledge, there are few studies jointly 
analyzing miRNAs and mRNA profiles to explore the 
roles of miRNAs on the cold adaptation of fish. Such a 
study would extend our understanding of miRNAs and 
environmental adaptation and undoubtedly lead to a more 
challenging and interesting question: whether and how 
miRNAs have integrated into cold-adaptive processes and 
contributed to the unique physiology of fishes endemic to 
extreme temperature conditions. 

Conclusions 

In this study, we systematically characterized the miR- 
NAome of zebrafish brain under normal and cold-accli- 
mated conditions using Solexa sequencing. Our data has 
greatly extended the zebrafish miRNAs repertory from the 
214 unique, mature miRNAs (miRBase databse, release 16 
version) to more than 574 ones. Besides, our results 
showed that the miRNAs, either individually or together, 
either directly or indirectly play only a minor role in tran- 
scriptomic plasticity during cold acclimation in the zebra- 
fish brain, suggesting that the cold-stress response of 
mRNA expression may be governed mainly through regu- 
latory modes other than post-transcriptional regulation by 
miRNAs. MiRNAs in animal brains may act more as devel- 
opmental regulators than thermal adaptability regulators. 

Methods 

Zebrafish treatment and brain RNA preparation 

Adult zebrafish (Danio rerio) were alternatively divided 
into two groups: a control group and an experimental 
group. In each group, the sex ratio was 1:1. The control 
group was raised under the normal laboratory breeding 
temperature (28°C), and the experimental group was sub- 
jected to cold (10°C) stress. An abrupt transfer of zebrafish 
from normal to cold (10°C) environments will cause most 
of them to die. To prevent this, the temperature of the 
cold group was decreased from 28°C to 10°C at a rate of 
one degree per hour. After 10 days of cold adaption, the 
brains of the zebrafish were dissected, and total RNAs 
were extracted from these brains using the total RNA iso- 
lation kit (Promega, Madison, USA). All experimental 
research on animals was conducted with the approval of 
the Animal Research Ethics Committee of our institute. 

Microarray profile 

To evaluate the genome-wide change of gene transcrip- 
tion in cold stress, we performed hybridization analysis 



using the Zebrafish Genome Array (Affymetrix) as pre- 
viously described [70]. Briefly, mRNAs were isolated 
from 10 \ig of total RNAs derived from pooled zebrafish 
brains of acclimated and control samples using the 
Eukaryotic Poly-A RNA Control Kit (Affymetrix) in an 
RNase-free environment. Subsequently, mRNAs were 
subjected to double-strand cDNA synthesis using a one- 
cycle cDNA synthesis kit (Affymetrix). The resulting 
double-strand cDNAs were in vitro transcribed into bio- 
tin-labeling cRNAs using the IVT Labeling kit (Affyme- 
trix). After that, biotin-labeling cRNAs were fragmented 
by heating and hybridized with the Affymetrix Zebrafish 
Genome Arrays. Finally, the signal intensity of the chip 
was scanned using the Gene Chip Scanner 3000 (Affy- 
metrix) and analyzed with the Gene Chip Operating 
Software (Affymetrix). The data derived from the gen- 
ome-wide expressional analysis was deposited in the 
GEO bank with the accession number GSE32092. 

Small-RNA deep sequencing 

For small-RNA Solexa sequencing, 20 (ig of total RNAs 
were subjected to electrophoresis on a polyacrylamide 
gel under denaturing conditions. The small RNAs ran- 
ging from 18-30 nt were excised and purified. The 
resulting small RNAs were sequentially ligated to a 5' 
adaptor and a 3' adaptor. The products were then 
reverse transcribed (RT) and amplified by polymerase 
chain reaction (PCR). The RT-PCR products of approxi- 
mately 90 bp (small RNA + adaptors) were isolated and 
subjected to sequencing analysis using the Illumina 
Genome Analyzer (Illumina, San Diego, USA) according 
to the manufacturer's instructions. 

MiRNA identification and annotation 

High-quality, small-RNA reads larger than 18 nt were 
extracted from raw reads and mapped to zebrafish gen- 
ome sequences (Danio_rerio.Zv8.56.dna.toplevel.fa.gz) 
using SOAP2 [37]. The miRNAs and their precursor 
structures were predicted by MIREAP (https://source- 
forge.net/projects/mireap/) under the default settings 
using the perfect-matched reads. Like the miRDeep soft- 
ware [39], MIREAP integrates miRNA biogenesis, 
sequencing depth and structural features to identify 
miRNAs and further interrogate their expression level 
from deep sequenced small RNA libraries. Stem-loop 
hairpins were retained only when they comply with: 1) 
the mature miRNAs-associated reads are mapped in the 
arm region of the precursors; 2) the free energy of the 
secondary structure calculated by RNAfold [71] is lower 
than -18 kcal/mol. In detail, the potential mature 
miRNA is defined as the most abundant read sequence 
that aligns to the potential precursor sequence; the 
expression level of a miRNA is then specified as the 
sum of an ensemble of reads that align with the 
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potential mature molecules, allowing three nucleotides 
sliding beyond the position of the potential mature 
miRNA at the 5' end. We then performed blast analyses 
on the predicted miRNAs precursors and mature miR- 
NAs against sequences deposited in the miRBase 
(release 16) database to annotate these miRNAs. Novel 
miRNAs were further tagged as "exon-derived", "intron- 
derived" or "intergenic" miRNAs according to the loci of 
their precursors on the genome. In addition, MiPred 
[41] was used to distinguish between genuine and 
pseudo miRNA precursors. This software incorporates 
the local, contiguous structure-sequence composition 
and the minimum free energy (MFE) of the secondary 
structure to eliminate potential false positive predictions. 

Frequency-based small-RNA reads clustering 

Following the method proposed by Qu et al. [43], the 
raw reads were clustered according to their read num- 
bers. After clustering, the accuracy of mapping short 
reads can be improved. 

Gene expression validation by q RT-PCR analysis 

Total RNA (8 |ig) from the cold and normal fishes was 
respectively reverse transcribed by Takara PrimeScript RT- 
PCR Kit (Takara, CN) with OligodT primers, according to 
the protocol the manufacturer recommended. Another 4 
(ig total RNA from the two treatments was used for small 
RNA (< 200 bp) isolation by Ambion mirVana miRNA Iso- 
lation Kit. Poly (A) tail was added to the small RNAs by 
Ambion Poly (A) Polymerase (Ambion, CN). 2 pairs of pri- 
mers for each selected gene or miRNA were designed and 
tested for amplification efficacy and specificity by agarose 
electrophoretic analysis of the PCR products. The primer 
pair with the best specificity was selected for further use in 
q-PCR. Q-PCR was run by using the SYBR RT-PCR kit in 
a Bio-rad CFX 96 (Bio-rad) machine and results were ana- 
lyzed by the CFX Manager software. The PCR thermal 
cycling was performed using 45 cycles of 95°C for 15 sec- 
onds, 72°C for 30 seconds and 60°C for 30 seconds. The 
PCR reaction for each gene was performed in triplet with 
either the housekeeping gene Beta-actin (for arrayed genes) 
or U6 RNA (for miRNAs) as control. 

Word enrichment/depletion in 3' UTR 

The mRNAs were sorted according to their significance 
of expression levels and their corresponding 3' UTRs. 
The longest sequences were retrieved from the Ensembl 
website using Biomart tools [49]. The word enrichment/ 
depletion patterns of these genes were generated by 
Sylamer [53]. 

Regulatory effect score 

The regulatory effect score (RE score) was defined as the 
inhibitory effect of a miRNA on the mRNAs in a 



sample, according to the average difference in expres- 
sion of its target versus non-target mRNAs [52]. For 
each miRNA, we first calculated the rank-based RE 
score for every sample and then compared the differ- 
ence between the normal and cold-acclimated groups 
using a two-sample t-test. The FDR values were calcu- 
lated using the R multtest [72] package. 

GO enrichment test 

For the list of mRNAs, we tested whether each had 
enriched GO terms in biological processes and molecu- 
lar functions using the GO:TermFinder package [46]. 
The GO annotation file was downloaded from http:// 
zfin.org/zf_info/downloads.html#go/on January 29, 2010. 
The GO ontology file was downloaded from http:// 
www.geneontology.org/ on January 23, 2008. 

miRNA target predictions 

Currently, there are several popular tools for predicting 
miRNA targets, including PITA [45], PicTar [73], miR- 
anda [51], mirWIP [74] and TargetScan [75]. These 
tools can be generally divided into two classes, methods 
that integrate the information of sequence conservation 
across multiple genomes and those that do not take this 
into account. Because of the limited genome data for 
zebrafish relatives, we predicted the miRNA targets by 
two methods, PITA and miRanda, independently. PITA 
scores the miRNA-target interactions based on a ther- 
modynamic model. A salient feature of this program is 
that it considers the site accessibility of the surrounding 
mRNA sequence but disregards sequence conservation 
between genomes. Thus, non-conserved or even species- 
specific target genes could be identified using this 
method. For each miRNA, the target genes are defined 
as those predicted to contain at least one binding site 
(a ag < -10) at their 3' UTR. We also ran miRanda 
(August 2010 Release) to predict the target genes with 
default parameters except that the energy parameter is 
set to -20. 

Network analysis 

The miRNA-mRNA relationship is modeled as a bipar- 
tite where an edge links a miRNA and an mRNA with a 
predicted binding site in its 3' UTR corresponding to 
the miRNA. We defined coherent and incoherent edges 
as those for which their associated miRNAs and mRNAs 
represented the opposite and parallel expression change 
under cold stress, respectively. We then calculated a 
regulatory density (RD) for a group of miRNAs and 
mRNAs within a network according to the following 
formula: 

RD = (#CE - #IE)/(#CE + #IE), (1) 
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where #CE (#IE) denotes the number of coherent 
(incoherent) edges. 

To estimate the significance of the observed RE, we 
generated 1,000 cohorts of random networks and recal- 
culated the RE by shuffling the original network while 
keeping the degree preservation according to the pub- 
lished method [76]. 

Additional material 
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